function outFeatures = genGrayFs(image)
%get 12 gray features
image = double(image);
image = image(:);
image(find(image<300)) = [];

grayBin = 100:10:3500;
grayHist = hist(image,grayBin);
%N(1:100) = 0;
%N(301:end) = 0;
grayHist = grayHist(100:300);
sumOfGrayHist = sum(grayHist);
grayHist = 100*grayHist./sumOfGrayHist;

grayTmp = grayHist;
grayTmp(find(grayTmp<0.25)) = [];
% grayTmp = grayTmp./sum(grayTmp);
outFeatures(1) = mean(grayTmp);
outFeatures(2) = std(grayTmp);
outFeatures(3) = skewness(grayTmp);
outFeatures(4) = Kurtosis(grayTmp);
sumx = 0;
sumxx = 0;

% Y(find(Y<1500)) = 0;
for j = 100:300
    sumx = sumx+grayHist(j-99)*(j-99);
    sumxx = sumxx + grayHist(j-99)*(j-99)*(j-99);
end
outFeatures(5) = sumx/sum(grayHist);
outFeatures(6) = sumxx/sum(grayHist) - outFeatures(5)*outFeatures(5);

%second 6 features
grayTmp2 = grayTmp.*grayTmp;
outFeatures(7) = mean(grayTmp2);
outFeatures(8) = std(grayTmp2);
outFeatures(9) = skewness(grayTmp2);
outFeatures(10) = Kurtosis(grayTmp2);
sumx = 0;
sumxx = 0;
%   Y(find(Y<1500)) = 0;
grayHist2 = grayHist.*grayHist;
for j = 100:300
    sumx = sumx+grayHist2(j-99)*(j-99);
    sumxx = sumxx + grayHist2(j-99)*(j-99)*(j-99);
end
outFeatures(11) = sumx/sum(grayHist2);
outFeatures(12) = sumxx/sum(grayHist2) - outFeatures(11)*outFeatures(11);